MODULE sefmm_cheb_multipoles

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: sefmm_build_cheb_multipoles,      &
             sefmm_contract_cheb_multipoles

   CHARACTER(LEN=10), PARAMETER :: Name = 'multipoles'
   INTEGER(INTK), SAVE :: IOwrite1
   INTEGER(INTK), SAVE :: IOwrite2
   REAL(REALK),   SAVE :: sefmm_multipole_threshold

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_build_cheb_multipoles(scheme,lmax,nbox,full_boxes)

      USE sefmm_density,  ONLY: sefmm_get_spectral_density
      USE fmm_car_to_sph, ONLY: fmm_transform_car_to_sph,      &
                                fmm_init_car_to_sph,           &
                                fmm_free_car_to_sph

      IMPLICIT NONE
      TYPE(scheme_paras), INTENT(IN)  :: scheme
      INTEGER(INTK),      INTENT(IN)  :: lmax, nbox
      INTEGER(INTK),      INTENT(OUT) :: full_boxes

      REAL(REALK), POINTER :: dens(:,:,:)
      REAL(REALK) :: Ine(1024+lmax, 0:lmax)
      REAL(REALK) :: car_mpoles((lmax+1)*(lmax+2)/2, 0:lmax)
      REAL(REALK) :: sph_mpoles(2*lmax+1, 0:lmax)
      REAL(REALK) :: T0, T1, grain

      INTEGER(INTK) :: i, npts, l,m

      grain = scheme%grain
      sefmm_multipole_threshold = scheme%sefmm_multipole_threshold

      CALL CPU_TIME(T0)
      WRITE(LUPRI,*) 'Writing file of boxed multipole moments...'

      CALL sefmm_open_mpoles_interface_files
      CALL sefmm_build_cheb_mpoles(lmax,grain,Ine)
      CALL fmm_init_car_to_sph(lmax)

      full_boxes = 0
      DO i = 1, nbox

         car_mpoles(:,:) = zero
         sph_mpoles(:,:) = zero
         CALL sefmm_get_spectral_density(i,dens,npts)
         IF (npts == 0) CYCLE
!print *, 'box =', i
!print *, npts, SIZE(dens,1), SIZE(dens,2), SIZE(dens,3)

         CALL sefmm_build_car_mpoles(lmax,npts,dens,Ine(1:npts,:),car_mpoles)

!      print *, '-------'
!      print *, 'cartesian'
!      DO l = 0, lmax
!         print *, '-------'
!      DO m = 1, SIZE(car_mpoles,1)
!         print *, l,m, car_mpoles(m,l)*half
!      end DO
!      end DO

         CALL fmm_transform_car_to_sph(car_mpoles,sph_mpoles,1,lmax)
!print *, '----------'
!print *, 'spherical'
!      DO l = 0, lmax
!print *, '----------'
!      DO m = 1, 2*l+1
!         print *, l,m, sph_mpoles(m,l)*half
!      end DO
!      end DO
!

         CALL sefmm_store_sph_mpoles(i,lmax,sph_mpoles,full_boxes)

      END DO

      CALL sefmm_close_mpoles_interface_files(lmax,full_boxes)
      CALL fmm_free_car_to_sph
      print *, "Number of boxes with multipole expansions =", full_boxes

      CALL CPU_TIME(T1)
      CALL TIMTXT('>>> TIME USED build cheb multipoles', T1-T0, LUPRI)

   END SUBROUTINE sefmm_build_cheb_multipoles

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_build_car_mpoles(lmax,npts,dens,Ine,car_mpoles)

      USE sefmm_cheb_utils, ONLY: cheb_get_Tn_ints

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: lmax, npts
      REAL(REALK),   INTENT(IN)  :: dens(0:npts-1,0:npts-1,0:npts-1)
      REAL(REALK),   INTENT(IN)  :: Ine(0:npts-1,0:lmax)
      REAL(REALK),   INTENT(OUT) :: car_mpoles((lmax+1)*(lmax+2)/2, 0:lmax)

      INTEGER(INTK) :: x,y,z, i,j,k, i0,j0,k0, l,m
      REAL(REALK)   :: tmp, fac

      DO l = 0, lmax
         m = 0
         DO x = lmax, 0, -1
            DO y = lmax, 0, -1
               DO z = lmax, 0, -1
                  IF ( x+y+z /= l ) CYCLE
                  m = m +1

                  ! Ine(n,e) = 0.0 for n+e = 2k+1
                  i0 = 1; IF ( MOD(x,2) == 0 ) i0 = 0
                  j0 = 1; IF ( MOD(y,2) == 0 ) j0 = 0
                  k0 = 1; IF ( MOD(z,2) == 0 ) k0 = 0

                  tmp = zero
                  DO k = k0, npts-1, 2
                     DO j = j0, npts-1, 2
                        fac = Ine(j,y)*Ine(k,z)
                        DO i = i0, npts-1, 2
                           tmp = tmp + dens(i,j,k)*Ine(i,x)*fac
                        END DO
                     END DO
                  END DO
                  car_mpoles(m,l) = tmp

               END DO
            END DO
         END DO
      END DO

   END SUBROUTINE sefmm_build_car_mpoles

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_build_cheb_mpoles(lmax,lambda,Ine)

      USE sefmm_cheb_utils, ONLY: cheb_get_Tn_ints

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: lmax
      REAL(REALK),   INTENT(IN)  :: lambda
      REAL(REALK),   INTENT(OUT) :: Ine(0:,0:)

      INTEGER(INTK) :: npts, e,n

      npts = SIZE(Ine,1)
      Ine(:,:) = zero

      ! Get I_{n}^{0} integrals to start recurrence
      CALL cheb_get_Tn_ints(npts-1,Ine(0:,0))

      ! Get I_{n}^{e} integrals by simple recurrence
      DO e = 1, lmax
         Ine(0,e) = Ine(1,e-1)
         DO n = 1, npts -1 -e
            Ine(n,e) = half*(Ine(n+1,e-1) + Ine(n-1,e-1))
         END DO
      END DO

      ! Rescale according to box size
      DO e = 0, lmax
         Ine(:,e) =  Ine(:,e) * (lambda*half)**(e+1)
      END DO

   END SUBROUTINE sefmm_build_cheb_mpoles

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_open_mpoles_interface_files

      IMPLICIT NONE
      CHARACTER(LEN=255) :: FBuf

      IOwrite1 = LUINTM +1
      IOwrite2 = LUINTM +2

      ! Open file of boxed-density multipole expansions
      FBuf = TRIM(Name)//".fmm1"
      OPEN(UNIT=IOwrite1, FILE=TRIM(FBuf), STATUS='REPLACE',    &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IOwrite1)

      ! Open file of box centres at which the FMM potential is computed
      FBuf = TRIM(Name)//".fmm2"
      OPEN(UNIT=IOwrite2, FILE=TRIM(FBuf), STATUS='REPLACE',  &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IOwrite2)

   END SUBROUTINE sefmm_open_mpoles_interface_files

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_close_mpoles_interface_files(lmax,nbox)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: lmax, nbox
      CHARACTER(LEN=255) :: FBuf

      ! Mark end of moments file with negative L-index
      WRITE(IOwrite1) 0,-1,0,0,0, 0d0,0d0,0d0, 0d0
      CLOSE(UNIT=IOwrite1, STATUS='KEEP')
      CLOSE(UNIT=IOwrite2, STATUS='KEEP')

      ! Write RHS boxed moments header file
      FBuf = TRIM(Name)//".fmm1header"
      OPEN(UNIT=LUINTM, FILE=TRIM(FBuf), STATUS='REPLACE',   &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      WRITE(LUINTM) lmax, 1, nbox
      CLOSE(UNIT=LUINTM, STATUS='KEEP')

      ! Write LHS potential grid points header file
      FBuf = TRIM(Name)//".fmm2header"
      OPEN(UNIT=LUINTM, FILE=TRIM(FBuf), STATUS='REPLACE',   &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      WRITE(LUINTM) nbox
      CLOSE(UNIT=LUINTM, STATUS='KEEP')

   END SUBROUTINE sefmm_close_mpoles_interface_files

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_store_sph_mpoles(box,lmax,mpoles,counter)

      USE fefmm_box_utils, ONLY: fefmm_box_centre

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)    :: box, lmax
      REAL(REALK),   INTENT(IN)    :: mpoles(2*lmax+1,0:lmax)
      INTEGER(INTK), INTENT(INOUT) :: counter

      REAL(REALK)   :: box_cntr(3)
      INTEGER(INTK) :: l,m, k
      LOGICAL :: flag

      box_cntr = fefmm_box_centre(box)

      flag = .FALSE.
      DO l = 0, lmax
         DO k = 1, 2*l +1
            m = k-l-1
            IF (ABS(mpoles(k,l)) > sefmm_multipole_threshold) THEN
               WRITE(IOwrite1) (counter+1),l,m,1,1, box_cntr(:), mpoles(k,l)
               flag = .TRUE.
            END IF
         END DO
      END DO

      IF (flag) THEN
         ! Grid points where multipole potential is computed
         WRITE(IOwrite2) 0,0,0,0,0, box_cntr(:), 1.0D0
         counter = counter +1
      END IF

   END SUBROUTINE sefmm_store_sph_mpoles

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_contract_cheb_multipoles(nbox,lmax,Vpot,energy)

      USE fmm_qlm_utils, ONLY: fmm_renormalise_qlm

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: nbox, lmax
      REAL(REALK),   INTENT(IN)  :: Vpot((1+lmax)**2,nbox)
      REAL(REALK),   INTENT(OUT) :: energy

      CHARACTER(LEN=255) :: FBuf

      REAL(REALK)   :: sph((1+lmax)**2,nbox)
      REAL(REALK)   :: box_cntr(3), tmp
      INTEGER(INTK) :: i, idum, l,m,lm, box

      ! File of boxed multipole moments
      FBuf = TRIM(Name)//".fmm1"
      OPEN(UNIT=LUINTM, FILE=TRIM(FBuf), STATUS='OLD',ACTION='READ',   &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(LUINTM)

      sph(:,:) = zero
      DO
         READ(LUINTM) box, l,m, idum,idum, box_cntr(:), tmp
         IF (l == -1) EXIT
         lm = l*(l+1) +m +1
         sph(lm,box) = tmp
      END DO
      CLOSE(UNIT=LUINTM, STATUS='KEEP')

      ! Renormalise moments before contraction with potential
      CALL fmm_renormalise_qlm(lmax,sph)

      energy = 0
      DO i = 1, nbox
         energy = energy + DOT_PRODUCT(Vpot(:,i),sph(:,i))
      END DO
      energy = half*energy

   END SUBROUTINE sefmm_contract_cheb_multipoles

!-------------------------------------------------------------------------------

END MODULE sefmm_cheb_multipoles
